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Abstract 

We report molecular simulation and experimental results for simple fluids 
adsorbed in activated carbon fibers (ACF), where the adsorbed phase con- 
sists of either one or two molecular layers. Our molecular simulations involve 
smooth pore-walls for large system sizes and a systematic system size-scaling 
analysis. We provide calculations of the Ginzburg parameter to monitor the 
self-consistency of our finite size simulation results, based on which we estab- 
lish that for system sizes smaller than 60 molecular diameters, fluctuations 
are too large to uphold the finite size simulation results. We present scaling 
analysis and free energy results for a system size of 180 molecular diameters. 
Our scaling analysis reveals a two-stage melting in our bilayer system con- 
sistent with KTHNY scaling for the order parameter correlation functions. 
Based on the Lee-Kosterlitz scaling of the free energy surface, we establish 
that the transitions are first order in the thermodynamic limit. Similar results 
for a monolayer established that the transitions were continuous. We provide 
scaling arguments to suggest that the change in the order of the transitions in 
going from a monolayer to bilayer adsorbed system results from the interac- 
tions of defects in the multilayers. We also report experimental measurements 



of CCI4 and aniline (with two confined molecular layers) adsorbed in ACF. 
The differential scanning calorimetry and dielectric relaxation spectroscopy 
measurements for the transition temperatures are in quantitative agreement 
with each other and with our simulation results. We also report nonlinear 
dielectric effect (NDE) measurements and show that our data is consistent 
with the NDE scaling law for the KTHNY scenario, which we derive. Fi- 
nally, we discuss the significance of the six-fold corrugated potential and the 
effect of the strength of pore-wall potential on the melting behavior. Both 
the simulations and experiments show that the hexatic phase is stabilized by 
confinement in carbon pores; for carbon tetrachloride and aniline the hexatic 
phase is stable over a temperature range of 55 K and 26 K, respectively. More- 
over, this stability increases as the ratio of fluid-wall to fluid-fluid attraction 
increases. 
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I. INTRODUCTION 



For continuous symmetry breaking transitions (such as melting transitions) in two- 
dimensions, the Mermin- Wagner theorem states that true long range order ceases to exist 1 . 
Halperin and Nelson proposed the "KTHNY" (Kosterlitz-Thouless-Halperin-Nelson- Young) 
mechanism for melting of a crystal in two dimensions 2 , which involves two transitions of 
the Kosterlitz-Thouless (KT) kind 3 . The crystal to hexatic transition occurs through the 
unbinding of dislocation pairs, and the hexatic to liquid transition involves the unbinding 
of disclination pairs. Each KT transition is accompanied by a non-universal peak in the 
specific heat above the transition temperature, associated with the entropy liberated by the 
unbinding of the vortex (dislocation or disclination) pairs, and by the disappearance of the 
stiffness coefficient associated with the presence of quasi-long-range order in the system. The 
KTHNY theory predicts that the correlation function associated with the translational order 
parameter in the crystal decays algebraically with exponent rj < 1/3, while long range orien- 
tational order is maintained, and the correlation function associated with the orientational 
order parameter in the hexatic phase decays algebraically with exponent < i]q < 1/4 while 
there is no translational order. The KTHNY theory is an analysis of the limit of stability of 
a two-dimensional solid, since it neglects the existence of the liquid phase and thereby does 
not impose an equality of the chemical potentials of the solid and liquid phases as the cri- 
terion for melting. Therefore, other pathways for two-dimensional melting can not be ruled 
out. For example, it is possible for the dislocation unbinding transition to be pre-empted 
by grain-boundary-induced melting, as shown by the work of Chui 4 , which predicts that 
the critical value of the defect core energy (E c ) above which the melting cross-over from 
grain boundary induced melting to two-stage KTHNY melting is E c = 2.8/c_bT 5 . Excellent 
reviews are available on the subject of two-dimensional melting 5 (also see section II). 

In this paper we report molecular simulation and experimental results for simple fluids 
adsorbed in activated carbon fibers (ACF), whose pore widths are such that they can ac- 
commodate one or two molecular layers. Our computer simulations involve smooth walls for 
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large system sizes, and we perform a systematic scaling analysis. We provide calculations 
of the Ginzburg parameter to monitor the self-consistency of our finite size simulation re- 
sults for system sizes (box lengths) up to 180 molecular diameters. We also present scaling 
analysis and free energy results for a system size of 180 molecular diameters. Based on the 
Lee-Kosterlitz scaling of the free energy surface, we establish the nature (first order vs. con- 
tinuous) of the transitions. We also report experimental measurements based on differential 
scanning calorimetry, dielectric relaxation spectroscopy, and nonlinear dielectric effect for 
CC1 4 and aniline (with two confined molecular layers) adsorbed in ACF. 

II. PREVIOUS WORK 

Experimental and computer simulation studies on the subject of two-dimensional melting 
often have the underlying objective of establishing whether or not the pathway of melting 
conforms to KTHNY behavior: (1) Does melting occur in two stages mediated by a hexatic 
phase? (2) Do the order parameter correlation functions associated with the crystal and 
hexatic phases have the appropriate scaling behavior? (3) Are the observed phase transitions 
first order or second order in the thermodynamic limit? The different experimental systems 
explored in this regard include free-standing liquid- crystalline (LC) films 6,7 , confined colloidal 
suspensions 8 ^ 11 , and adsorbed fluid on a planar substrate 12 ^ 19 . 

Computer Simulation Studies of Two-Dimensional Melting. Early simulation studies on 
small (< 10000 atoms), strictly two-dimensional systems failed to provide compelling evi- 
dence to support the KTHNY melting scenario 5,20-22 . For systems with repulsive interac- 
tions, it was shown by Bagchi et. al. 23 , using a systematic scaling analysis on large system 
sizes (~ 64,000 atoms), and subsequently by Jaster 24 , that the equilibrium properties are 
indeed consistent with the KTHNY theory of melting (a single stage melting via first-order 
transitions with finite correlations was ruled out). Bagchi et. al. 23 have also observed a 
first order single stage melting for a two-dimensional system of disks interacting with a soft 
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repulsive potential for a small system size, crossing over to two-stage continuous melting 
consistent with KTHNY scaling for larger system sizes, suggesting the earlier studies were 
plagued by serious system size effects. 

Computer simulation studies that directly mimic the liquid-crystalline thin film experi- 
ments have not been attempted owing to the complexity of the molecular architecture and 
experimental conditions. Idealized simulations of spheres confined between hard walls and 
interacting with each other via a model potential that mimic the screened potential in col- 
loidal systems 9 (but otherwise neglect the solvent) have been performed by Bladon and 
Frenkel 25 , and Zangi and Rice 26 . The former study reported a strong dependence (and 
qualitatively different phase diagrams of the melting region) of melting behavior on the 
parameters of the inter-molecular interaction of a two-dimensional square-well fluid, and 
attributed the difference to the changing values of the defect core energy 25 ; in particular, 
the study puts the two-dimensional hard disks in the single-stage first order melting regime 
(rather than the KTHNY regime). In subsequent experiments and simulations carried out by 
Rice and co-workers using two different screened (effective) potentials with which confined 
colloidal spheres interact 9 ' 10 ' 26 , a strong dependence of the melting scenario on the inter- 
particle interaction, consistent with the computer simulation study by Bladon and Frenkel 25 
was reported. A clear picture regarding the differing scenarios has not yet emerged due to: 
(a) the experimental studies on colloidal systems are faced with the question of attainment 
of thermodynamic equilibrium 9,10 , while the simulation results are for a small system size 25 , 
and are therefore subject to finite-size effects; (b) the calculated value for the defect core 
energy E c for hard disks lie in the KTHNY regime 27 ; (c) the qualitative behavior of these 
studies in the limit of a hard sphere potential differs from that of a recent simulation study 
by Jaster 24 on a large 2-d system of hard disks. 

Several computer simulations aimed to mimic an adsorbed fluid on graphite have been 
reported 28-30 . Abraham 28 has reported several studies of xenon on graphite in which the 
melting transition temperatures are in quantitative agreement with experiment. However, 
the existence of an intrinsic hexatic phase could not be demonstrated by the simulations, 
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which were performed on fairly small system sizes. Similar conclusions were drawn by Roth 29 
based on simulations of krypton on graphite, again for a small system size. Alavi 30 reported a 
simulation study of CD 4 on MgO (using a small system size) with melting behavior consistent 
with the KTHNY scenario and a defect core energy of 9ksT (KTHNY regime). One notable 
difference between the CD 4 on MgO system and the Kr, Xe on graphite systems is that, 
in the former case, due to the molecular structure of CD 4 adsorbed on MgO, the relevant 
degrees of freedom at the temperature of study make the adsorbate move in a field-free 
manner. 

In summary, the simulations on small systems have failed to provide compelling evidence 
of the melting picture, and therefore can not be used conclusively to fill in the gaps or 
interpret experimental measurements. While the system-size scaling analysis for large system 
sizes in the simulations of Bagchi et. al. 23 (which indicated a KTHNY melting scenario) 
clearly underlines the need to use large system sizes in studying the two-dimensional melting 
problem, even such methods do not provide a source of distinction between the first-order 
and continuous nature of the transition. 

Two-dimensional Melting in Adsorbates Confined in Porous Media. Activated carbon 
fibers (ACF) possess micro-crystallites made up of graphene sheets that tend to align in 
similar directions, with slit shaped voids between the microcrystals. The spontaneous order- 
ing of the molecules adsorbed in these voids into distinct two- dimensional molecular layers 
(analogous to the structure of a smectic-A phase in liquid crystals) makes the adsorbed phase 
a quasi-two-dimensional system. Microporous ACF can be prepared having a range of mean 
pore sizes, ranging from those that can accommodate just a single layer of the adsorbate in 
the micropores to those that can accommodate a few layers 31 . The microcrystals in ACF are 
themselves arranged as an amorphous matrix; therefore the maximum in-plane correlation 
length of the adsorbed fluid is thought to be limited by the average size of the microcrystals. 
Real samples of ACF have different average sizes for their microcystals (ranging from 5 nm 
to 100 nm) depending on the method of activation and the source of the carbon material 
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used to derive the porous material 31 . Electron micrographs are commonly employed to de- 
termine the average size of the microcrystals, while nitrogen adsorption isotherms are used 
to determine the average pore size and pore size distribution. The substrate field induced by 
the porous matrix is thought to have a six-fold symmetry like that of graphite. Additionally 
due to the high density of carbon atoms in the graphene microcrystals, the confined fluid 
feels a large potential well, and this increases in depth for decreasing pore sizes; this unique 
feature induces pronounced layering effects even in multi-layer adsorbate phases, each layer 
effectively being two-dimensional. The quasi-two-dimensional phase of adsorbed simple flu- 
ids (spherical molecules) do not have the complication of additional coupling due to any 
herring bone symmetry, in contrast to the LC systems. In addition, owing to the nanoscopic 
length scales, these systems are thermally driven and can be studied under equilibrium (a 
challenging issue for the colloidal systems). The large value of the defect core energy 30 in 
adsorbed systems and the role of out-of-plane motions in stabilizing the hexatic phase in 
these quasi-two-dimensional systems have been noted earlier 5,26 . 



III. MOLECULAR SIMULATION METHODS 

We performed Grand Canonical Monte Carlo (GCMC) simulations of a fluid adsorbed in 
slit-shaped pores of width H, where H is defined as the perpendicular distance between the 
planes passing through the nuclei of the first layer of molecules that make up the pore walls 
of the slit-shaped pore. The interaction between the adsorbed fluid molecules is modeled 
using the Lennard- Jones (12,6) potential with size and energy parameters, <jg, ande ff . The 
Lennard - Jones potential was cut-off at a distance of 5crg, beyond which it was assumed to 
be zero. The pore walls were modeled as a continuum of LJ molecules using the "10-4-3" 
Steele potential 32 , given by, 

m*> - aw** [i (^r - (^) 4 - ( 3A(z ; L 0MAr )} (i) 

Here, the cr's and e's are the size and energy parameters in the Lennard- Jones (LJ) potential, 
the subscripts f and w denote fluid and wall respectively, A is the distance between two 



successive lattice planes of pore wall, z is the coordinate perpendicular to the pore walls and 
p w is the number density of the wall atoms. For a given pore width, H, the total potential 
energy from both walls is given by, 

0porc(^) = <frfw(z) + 0f w (.ff - z) (2) 

The strength of attraction of the pore walls relative to the fluid-fluid interaction is determined 
by the coefficient a, 

a = P- e ^ A (3) 

in equation 1. The fluid- fluid interaction parameters were chosen to be those for CCI4 33 , with 
erg = 0.514 nm, e^/ks = 395 K; with these parameters the freezing temperature of bulk CCI4 
is reproduced. The parameters for the graphite pore wall interaction to mimic the micropores 
in ACF were taken from Steele 32 , (o ww = 0.34 nm, p w = 114 nm~ 3 , e ww //c£ = 28 K, and 
A = 0.335 nm). The Lorentz-Berthlot mixing rules, together with the ff and ww parameters, 
were used to determine the values of o fw and e fw . The simulation runs were performed by 
fixing the chemical potential, //, the volume ,V, of the pore and the temperature, T. Two 
pore-widths H = lo// + o/ w = 0.91 nm (monolayer of adsorbate) and H = 2<7ff + <jf w = 1.41 
nm (bilayer of adsorbate) , were chosen for study. The rectilinear dimension of the cells were 
therefore 180off x 180off x H (93 nm x 93 nm x H nm). Typically the system consisted 
of up to 64,000 adsorbed fluid molecules. Periodic boundary conditions were employed in 
the x and y directions (the xy-plane being parallel to the pore walls) and no long range 
corrections were applied. The adsorbed molecules formed distinct molecular layers parallel 
to the plane of the pore walls. The simulation was set up such that insertion, deletion and 
displacement moves were chosen at random with equal probability. The calculations for the 
larger of the chosen pore-width also enable a direct comparison to be made with experimental 
measurements for CC1 4 confined in porous pitched-based activated carbon fiber ACF A-10, 
of mean pore width H = 1.4 nm. We expect the approximation of a structureless graphite 
wall to be a good one here, since the diameter of the LJ molecule (0.514 nm) is much larger 
than the C-C bond length in graphite (0.14 nm). The state conditions were such that the 
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confined phase was in equilibrium with bulk LJ CCI4 at 1 atm. pressure. The simulations 
were started from a well equilibrated confined liquid phase at T = 400 K, and in successive 
simulation runs the temperature was reduced. Equilibration was for a minimum of 11 billion 
steps; the standard deviation of block averages of total energy of the system was 10~ 4 and 
the rate of insertion was equal to that of deletion to a factor of 10~ 6 . 

Free Energy Determination. Motivated by the work of Frenkel and co-workers 34-36 , we 
employ the Landau free energy approach that was successful in our earlier studies 33 ' 35-37 . 
The Landau-Ginzburg formalism 38 involves choosing a spatially varying order parameter 
$(f), that is sensitive to the degree of order in the system. We use a two- dimensional bond 
orientational order parameter to characterize the orientational order in each of the molecular 
layers that is defined as follows 39 : 



where iVj, is the number of nearest-neighbor bonds and 9k is the bond angle (see below); 
the summation is over the (imaginary) bonds joining the central molecule to its nearest 
neighbors. ^ 6 j(p) measures the hexagonal bond order at position p in the xy plane within 
each layer j, and is calculated as follows. Nearest neighbors were identified as those particles 
that were less than a cutoff distance r nn away from a given particle, and belonged to the 
same layer. We used a cutoff distance r nn = 1.3 as, corresponding to the first minimum of 
g(r). The orientation of the nearest neighbor bond is measured by the 9 coordinate, which is 
the angle that the projection of the nearest-neighbor vector on to the xy-pl&ne makes with 
the x axis. ^/ 6 j(p), is calculated using equation 4, where the index k runs over the total 
number of nearest neighbor bonds N b at position p , in layer j. The order parameter \I> 6 j 
in layer j is given by \I/ 6 j = | J dp^ % j (p ) | / J dp. For the case of L J CC1 4 in slit-shaped 
pores, where there is significant ordering into distinct molecular layers, the spatially varying 




(4) 
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order parameter $(r ) can be reduced to &(z), and can be represented by, 

n 

<S>(z) = £ * 6J <^-%) (5) 

In equation (5), the sum is over the number of adsorbed molecular layers, and % is the 
z coordinate of the plane in which the coordinates of the center of mass of the adsorbed 
molecules in layer j are most likely to lie. It must be recognized that each of the ^6j's are 
variables that can take values in the range [0, 1]; the number of layers is n = 1 or n = 2 in our 
system. The histograms are collected to evaluate the probability distribution P[^ 6:1 , ^6,2] ■ 
The Landau free energy 38 A[\l/ 6i i, \l/ 6i2 ] (for the bilayer system) is given by, 

A[* 6 ,i,* 6)2 ] = -k B T ln(P[* 6jl ,* 6)2 ]) + Constant (6) 

The Landau free energy is computed by a histogram method combined with umbrella 
sampling 40 , using the probability distribution P[^ 61 , ^6,2]- Detailed procedures to collect 
statistics, construct the histograms and to choose the the weighting functions for performing 
the umbrella sampling are described elsewhere 34 36 . The grand free energy of a particular 
phase A, Qa = —ksT ln(S) (where S is the partition function in the grand canonical en- 
semble), is related to the Landau free energy by, 

2 r 

exp(-/3fi) = J} / d^ 6 ,j exp(-/3A[tt 6il ,tt 6)2 ]) (7) 
3=1 J i 

where the limits of integration in equation 7 are from the minimum value of (^6,1, ^6,2) to 
the maximum value of (^6,1, ^6,2) that characterize the phase A. The grand free energy is 
computed via numerical integration of equation 7. The two- dimensional, in-plane positional 
and orientational correlation functions (gj(r) and G§j{r) of layer j), were monitored to 
keep track of the nature of the confined phase. The positional pair correlation function 
is the familiar radial distribution function. The orientational pair correlation function is 
Gejip) = (^lj(0)^Qj(p)), where p — xe x + ye y . The above equations correspond to the 
two-layer case. For the one layer case, the Landau free energy reduces to A[^ 6i i] and the 
index j in equation 7 assumes a value of 1. 
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IV. EXPERIMENTAL METHODS 



Sample Preparation. The liquid samples were reagent grade chemicals, and were distilled 
twice prior to use in the experiment. The conductivities of the purified dipolar fluid samples 
were found to be less than 10~ 10 ohm _1 m _1 . The microporous activated carbon fiber (ACF) 
samples used were commercially available from Osaka gas company, Japan, with a pore 
size distribution of about 5% around the mean pore diameter 41 . ACF samples with an 
average pore width of 1.41 nm were used. The pore samples were previously characterized 
by Oshida and Endo by obtaining electron micrographs 31 ' 42 , and by Kaneko and co-workers 
using nitrogen adsorption measurements 41 . The characterization results for ACF showed 
that these amorphous materials consisted of uniform pores formed by graphitic microcrystals, 
with an average microcrystal size of 7-10 nm 31 . 

Differential Scanning Calorimetry (DSC). A Perkin-Elmer DSC7 differential scanning 
calorimeter was used to determine the melting temperatures and latent heats of fusion, by 
measuring the heat released in the melting of aniline and CCI4. The temperature scale 
of the DSC machine was calibrated using the melting temperature of pure fluid from the 
literature. The temperature scanning rates used for the melting and freezing runs varied 
from O.lK/min to 0.5K/min. The background of each raw DSC spectrum was subtracted, 
based on a second-order polynomial fit to the measured heat flow away from the signals 
of interest. The melting temperatures in the bulk and confined systems were determined 
from the position of the peaks of the heat flow signals, and the latent heats were determined 
based on the scaled area under these signals. The melting temperature was reproducible to 
within 0.5 °C for fluid adsorbed in pores. The latent heats were reproducible to within 5%. 

Dielectric Relaxation Spectroscopy (DRS). The relative permittivity of a medium, 
k* = K r — iKi, is in general a complex quantity whose real part n r (also known as the di- 
electric constant) is associated with the increase in capacitance due to the introduction of 



11 



the dielectric. The imaginary component Ki is associated with mechanisms that contribute 
to energy dissipation in the system, due to viscous damping of the rotational motion of the 
dipolar molecules in alternating fields. The latter effect is frequency dependent. The exper- 
imental setup consisted of a parallel plate capacitor of empty capacitance C Q = 4.2 pF. The 
capacitance, C, and the tangent loss, tan(<5), of the capacitor filled with the fluid between 
the plates were measured using a Solartron 1260 gain impedance analyzer, in the frequency 
range 1 Hz - 10 MHz, for various temperatures. For the case of the adsorbate in ACF, the 
sample was introduced between the capacitor plates as a suspension of ground ACF parti- 
cles of 0.1 mm mesh ACF particles in pure fluid. Due to the large conductivity of ACF, 
the electrodes were blocked by teflon. The relative permittivity is related to the measured 
quantities by: 

ts>r — ^ i ™% — \°) 

In equation (8), C is the capacitance, C Q is the capacitance without the dielectric and S is the 
angle by which current leads the voltage. The complex dielectric permittivity, 
is measured as a function of temperature and frequency. 

For an isolated dipole rotating under an oscillating electric field in a viscous medium, the 
Debye dispersion relation is derived using classical mechanics 43 , 

K = K oo,r + : ; : (9) 

1 + IUJT 

Here u is the frequency of the applied potential and r is the orientational (rotational) 
relaxation time of a dipolar molecule. The subscript s refers to static permittivity (low 
frequency limit, when the dipoles have sufficient time to be in phase with the applied field). 
The subscript oo refers to the optical permittivity (high frequency limit) and is a measure of 
the induced component of the permittivity. Further details of the experimental methods are 
described elsewhere 44 ' 45 . The dielectric relaxation time (molecular orientational relaxation 
time for dipolar molecules) was calculated by fitting the dispersion spectrum of the complex 
permittivity near resonance to the Debye model of orientational relaxation. 
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Nonlinear Dielectric Effect (NDE). The nonlinear dielectric effect (NDE) is denned as 
the permittivity change, Ae, of the medium in a strong electric field, E: 

NDE = = (10) 

where e# is the permittivity of the medium in a field E. The permittivity e is related to the 
dielectric constant K r by the relationship K r = e/e Q , e Q being the permittivity in vacuum. 
Note that within linear response (for weak fields), the permittivity is independent of the 
applied electric field; However, for strong fields, the most general form for the permittivity 
(consistent with field reversal invariance) can be written as €e = €e=o + bE 2 + cE A + • • • . 
By definition, NDE in equation 10 represents the first order (non-linear) response consistent 
with the general equation. The NDE was measured using the pulse method using rectangular 
millisecond pulses of the electric field with amplitudes ranging from 4 x 10 7 to 9 x 10 7 V/m. 
The separation between the invar electrodes in the measuring condenser was 2 x 10~ 4 m, and 
the changes in the capacitance were measured to an accuracy of 5 x 10~ 4 pF 46 . 

V. MOLECULAR SIMULATION RESULTS 

Our simulations were made to mimic two experimental systems of CC1 4 confined in ACF 
with one as well as two confined molecular layers. The orientational correlation function 
G 6 j(r) for equilibrated liquid, hexatic and crystal phases are shown in figure 1 . The long- Fig. 1 
range orientational order in the crystal phase, the algebraic decay of orientational order in 
the hexatic phase and the exponential decay of orientational order in the liquid phase are 
captured in the plots. Since the phases are characterized by long-ranged correlations, we 
check for the attainment of equilibrium and artifacts due to finite size of the simulations as 
described below. 

Simulation results are always for a finite system size, and system-size effects can not be 
avoided. When the state conditions are such that the relevant order-parameter correlation 
length £ approaches or exceeds the spatial extent of the simulation cell L, the system crosses 
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over to a mean-field behavior 47,48 . For correlations in the bond-orientational order param- 
eter, this is the case in the hexatic phase at all temperatures, and in the liquid phase near 
the liquid-hexatic transition. Under these circumstances, the self consistency of the mean 
field result is checked by calculating the Ginzburg parameter 7gl 49 , 



(%j) 2 its,,) 2 y ' 

where L 2 x H = V, the volume of the system, L is the length of the simulation box, and 
I(L) is given by 

I(L) = k (12) 

Jv d P 

The Ginzburg parameter gives the ratio of the variance in the order parameter (due to 
thermal fluctuations) to the square of the average value of the order parameter. Mean field 
results are clearly suspect in the regime 7ql ~ 1, when long-range fluctuations are likely 
to destroy the observed order in the finite size simulations. If 7 G l <^ 1, the simulations are 
self-consistent, and the ordered phase observed in the simulations is stable against thermal 
fluctuations. The Ginzburg parameter for different system sizes is calculated (using equa- 



tion 11) by numerically integrating equation 12, and are reported in Table I. Our calculations | Tab I 
for the bilayer system (H = 1.41 nm) indicate that the simulation results for these quasi- 
two-dimensional systems are only reliable for systems sizes > 60er ff . For smaller system 
sizes, the bond-orientational order parameter fluctuations are too large to validate the finite 
size results from computer simulations. Therefore, we report the free energy calculations 
and scaling behavior for system sizes of 60 and 180<r ff . 

To quantify the scaling properties of the orientational correlation function (figure 1) in 



each phase, we plot log [/(Lb) /J(L)] as a function of log [Lb/L] in figure 2, for different Fig. 2 
values of the block length. The block length L B < L defines the subsystem L B x L B x H for 
which the integrals (equation 12) were evaluated; the integrals I {Lb) and II were calculated 
by numerically integrating equation 12. The expected KTHNY scaling for liquid, hexatic, 
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and crystal are given by 2 : 

exp(— r/£) liquid 
G ej( r ) ~ \ r- 1 ' hexatic ( 13 ) 

Constant 7^ crystal 

where £ is the correlation length and 77 is the exponent characterizing the algebraic decay of 
the order parameter correlation function. (The exponent 77 has its origins from the Ornstein- 
Zernikie formalism 50 discussed later, however, the exponent values and functional forms for 
liquid-liquid critical phenomena and KTHNY formalisms are different from one-another) . 
Equations 12 & 13 lead to the following scaling for I(L): 

Constant / L 2 liquid 
I(L) = { l~7(2 - n) hexatic ( 14 ) 
Constant 7^ crystal 

Therefore, in figure 2, a slope of —2 corresponds to the liquid phase having a finite correlation 
length; a slope of —1/4 corresponds to an algebraic decay of G 6j (r) with exponent —1/4 
and is observed for the hexatic phase near the hexatic-crystal transition; a slope of is 
observed for the crystal phase, confirming the true long-range orientational order. The 
scaling properties shown in figure 2 is further evidence for the attainment of equilibrium in 
our simulations. 

The order of the phase transition was determined by examining the dependence of the 
free energy barrier separating two phases as a function of system size. In doing so, we 
consider only those system sizes for which 7ql « 1 In order to determine the order of 
the phase transition using the Lee-Kosterlitz scaling analysis 51 of the free energy surface, 
we obtained the Landau free energy functions at the exact transition temperatures for the 
corresponding system size. This was done as follows: the Landau free energy surface was 
calculated at a temperature close to the transition, from which we calculated the grand 
free energies at that temperature. Then, by numerically integrating the Clausius-Claperon 
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equation, d(Q/T)/d(l/T) = (U) — n(N), we located the exact transition temperature. The 
Landau free energy function was re-calculated at the transition temperature by using a 
weighting function equal to exp(/3A[\l/ 6i:; ]) from the initial calculation. This procedure was 
repeated for each system size. The Landau free energy and the grand free energy results 



for the system size of 180(7^ are provided in figure 3. The relative stability of the liquid, 
hexatic and crystalline phases are inferred from the Landau free energy plot (figure 3), and 
the exact transition temperatures are given by the grand free energy plot (figure 4). For the 



Fig. 


3 




Fig. 


4 



bilayer adsorbed fluid, the free energy barrier separating hexatic and crystal phases is a linear 
function of system size (figure 3a). This implies that the mechanism of phase transition is 
via nucleation and is a clear indication of a first order phase transition in the thermodynamic 
limit. The critical nucleus at the transition temperature is equal to the system size, implying 
a L d_1 (d being the dimensionality of the system) dependence on system size 51 . In contrast 
to the bilayer case, the Lee-Kosterlitz scaling of the free energy surface for the adsorbed 
monolayer (H = 0.91 nm) in figure 3b establishes that the transitions are continuous in the 
monolayer. In this case the free energy barrier separating the liquid and the hexatic phase 
is independent of system size, a clear signature of second order transition 51 . We attribute 
this remarkable change in the mechanism of phase transition (spinodal decomposition in the 
monolayer to nucleation in the bilayer) to the interaction of defect structures across layers 
(see section VII). 



VI. EXPERIMENTAL RESULTS 



The DSC and DRS results for CC1 4 and aniline confined in ACF are provided in fig- 



ures 5 and 6. The DSC scans show two peaks for each fluid reminiscent of liquid-hexatic Figs. 5,6 
(high temperature) and hexatic-crystal (low temperature) heat capacity peaks associated 
with Kosterlitz-Thouless phase transitions. The DRS measurements for the dielectric con- 
stant for confined aniline confirm the existence of two phase transitions (figure 6a). In 
addition, the molecular orientational relaxation times r, (figure 6b) are consistent with the 
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existence of liquid (r ~ns), hexatic (r ~ /is), and crystalline (r ~ms) phases. The transition 
temperatures from simulations and DSC for CC1 4 are in good agreement; also, the transition 
temperatures inferred from the DSC and DRS results for aniline are in agreement (Table II). 



The NDE measurements for confined CC1 4 and aniline in ACF are provided in figure 7 and Fig. 7 
Table II. The NDE signals show signatures of divergence at the liquid-hexatic and hexatic- Tab II 
crystal transition temperatures. The NDE results for the transition temperatures for CCI4 
and aniline are in near quantitative agreement with the simulation result, and also with those 
from DSC, and DRS (Table II). Below, we show that the scaling of the NDE signal with 
temperature is consistent with the KTHNY theory for liquid-hexatic and hexatic-crystal 
transitions. 

It was empirically known that the dipolar fluctuations associated with polarizability con- 
tributed to a positive signal for NDE. The anomalous increase in NDE in the vicinity of the 
phase transition point, reflecting the onset of long-range correlations and associated dipo- 
lar fluctuations in the system, has been experimentally documented 45 . Based on a droplet 
model for treating critical fluctuations by Oxtoby and Metiu 53 , the dipolar fluctuations near 
a critical point have been correlated (linearly) with the relevant order parameter fluctua- 
tions. De Gennes 54 showed that the conformity to the droplet model by Oxtoby and Metiu 53 
immediately leads to NDE scaling laws that are similar to those for critical fluctuations in 
light scattering experiments. The scaling law for critical opalescence in light scattering 
(where the order parameter correlation function is g(r) associated with density fluctuations) 
is given by the familiar expression 49 , 



N r 

m = < £ e M ii ■ exp W - x$)) =N 1 + y 9 (r) exp( if . f) #r. (15) 



a,0=l 

In the above equation, I(q) = NS(q), is the intensity measured in scattering experiments, 
S(q) is the familiar structure factor, d is the number of dimensions and N the total number 
of molecules. The double summation is over the total number of molecules, where x a ,X/3 
represent the coordinate vector of the given molecule. The second equality in equation 15 fol- 
lows from evaluating the ensemble average 49 . By employing the extended Ornstein-Zernikie 
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representation 50 for the density-density correlations (above T c ), namely, 

9(r) ~ 1 ~ exp(-r/0, (16) 

where £ is the correlation length. The integral in equation 15 in the limit q — > reduces to 

HQ) ~ i d ~ l ~ V - (17) 

The equation 16 corresponds to Fischer's extension of the original Ornstein-Zernikie formal- 
ism to alleviate a logarithimic divergence of g(r) for large r in two-dimensions. Therefore, 
the equation provides the definition for the critical exponent r\. The exponent satisfies the 
equality {2—rf)v = 7, for systems with Ising symmetry, where v (defined below) is the critical 
exponent associated with the correlation length, and 7 that with isothermal compressibility. 
According to the droplet model 53 , the analogous expression for NDE is given by, 

NDE ~ J ($(0)$(f))exp(iq-f)d d r, (18) 

where, (<&(0)$(f)) is the order parameter correlation function and d the dimensionality of 
the system. The rationale behind the droplet model is that the NDE signal results from the 
scattered intensity caused by the effective dipole-dipole correlations in the system, the same 
way light scattering intensity results from density-density correlations. 

As shown by de Gennes and Prost 54 , the droplet model leads to the experimentally 
observed scaling for NDE in the vicinity of a liquid-liquid critical point in three dimensions, 
namely, NDE ~ f^ 1 ^', with f = £ t~ u , and t = \T - T c \/T c . The theoretically predicted 
NDE scaling for this system (Ising symmetry, d=3) is therefore, 

d — 3, Ising; e.g., liquid-liquid 
v = 0.64, i] = 0.041; 3-d Ising model 
v = 0.64, 77 = 0; Ornstein-Zernikie ( i9 ) 
v = 0.5, rj = 0; Mean field theory 
dipolar fluids; experiment 45 . 
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Encouraged by the agreement of the phenomenological NDE scaling laws with experi- 
ment for systems with Ising symmetry, we derive the analogous NDE scaling laws near the 
liquid-hexatic, and hexatic-crystal phase boundaries. Our analysis was greatly facilitated 
because Halperin and Nelson 2 and Young 2 have worked out the scaling of the relevant order 
parameter correlations to be used in equations 15 and 17. The order parameter correlations 
above the transition temperatures for the liquid-hexatic and hexatic-crystal transitions have 
the form, Gt,q{j-) ~ exp(— r/£+), where the subscripts T, 6 correspond to translational and 
orientational correlation functions, and + refers to the correlation function above the tran- 
sition temperature, i.e., t > 0. Above the hexatic-liquid phase transition temperature, 
£ + = A exp(B/\t\ u ), with v = 0.5 based on a renormalization group analysis. Similarly, 
the translational order parameter correlation length is given by £ + = Aexp(B/\t\ u ), with 
v = 0.37 for a smooth substrate, above the crystal-hexatic phase boundary. The phenomeno- 
logical NDE scaling directly follows from using equation 17 (note that rj = 0, d = 2), 

1/NDE = Aexp(-B/|T-T c n (20) 

with v = 0.5 for the liquid-hexatic (high temperature) transition and v = 0.37 for the 
hexatic-crystal (low temperature) transition. 



The NDE scaling with temperature for confined aniline is re-plotted in figure 8 along Fig. 
with the NDE scaling law (equation 20). The results show the consistency of the NDE 
signals with KTHNY scaling laws for the liquid-hexatic and hexatic-crystal transitions. The 
estimates for the exponent v for liquid-hexatic and hexatic-crystal transitions based on the 
experimental measurements is currently not possible owing to the nature of the ACF sample; 
true divergence is not observed in our measurements as the length of the in-plane correlations 
in the fluid are curtailed by the finite size of the graphitic microcrystals in ACF. Therefore, 
NDE signals close to the transition temperatures for an infinite system are not accessible in 
our experiments. Instead, the NDE signals in our experiment correspond to a finite system 
size of 5 nm, the average size of the graphitic microcrystal (in the xy plane) in our ACF-A10 
sample. 
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VII. DISCUSSION 



Our computer simulations for a bilayer of adsorbed fluid in a slit-pore with smooth walls 
has established the existence of two-stage melting in the thermodynamic limit via system 
size scaling analysis. Based on the calculations of the Ginzburg parameter we establish that 
for system sizes smaller than 60 molecular diameters, fluctuations are too large to uphold the 
finite size simulation results. Based on the Lee-Kosterlitz scaling of the free energy surface, 
we establish that the transitions are first order in the thermodynamic limit. Similar results 
for a monolayer established that the transitions were continuous. Our experimental results 
rely on indirect evidence based on phase transitions (structural measurements are considered 
direct evidence). Nevertheless, the quantitative agreement between the simulations and 
DSC, DRS, and NDE measurements, make a strong case for the existence of KTHNY melting 
and therefore a hexatic phase in the confined system. This interpretation, if correct, shows 
that the hexatic phase is stable over a wide temperature range, 55 K for CC1 4 and 26 K for 
aniline. This large range of state conditions over which the hexatic phase is stable may be 
unique to confined fluids in porous media. 

Effect of the Six-Fold Potential. The main difference between the simulation results and 
the experimental results reported in this paper is that the former corresponds to a smooth 
substrate, while the latter corresponds to a substrate with 6-fold symmetric potential. On 
a theoretical basis, the effects of the six-fold symmetry of the pore potential and that of 
the finite size of the graphitic microcrystals combined might be expected to nullify the 
liquid-hexatic phase transition and cause a rounding of divergences associated with the 
hexatic-crystal transitions. The latter is observed in the NDE measurements (the NDE 
signals do not diverge), while the former manifests itself as a remnant KTHNY behavior, 
presumably due to a weak six-fold substrate potential. It is worth pointing out that the LJ 
diameters of CC1 4 and aniline (0.5-0.6 nm) are much larger than the C-C bond length in 
graphite (0.14 nm), so that the fluid molecules only feel a mild corrugation in the fluid- wall 
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potential in passing along the surface. This feature contrasts with the situation for earlier 
experiments concerning Xe, Ar, and Kr on graphite substrates, where the LJ diameters 
for the fluid are much smaller and comparable to the C-C spacing in graphite. There is 
evidence of an intrinsic hexatic phase even for the case of Xe in graphite 15 ' 16 ; therefore, 
for substrate molecules larger than Xe, we expect the crystal to melt into a hexatic phase 
with intrinsic stiffness with respect to bond-orientational fluctuations. While we note the 
quantitative agreement in transition temperatures between simulations and experiments 
and have described the respective domains of validity, the question of whether an intrinsic 
hexatic phase is observed in a realistic simulation with six-fold symmetric substrate potential 
remains to be answered. Simulations are currently in progress to address this question. On 
the experimental side, the challenge is to synthesize a porous sample with uniform pores 
with large microcrystal sizes to examine the scaling for several decades in \T — T c \/T c on 
a log scale, and to obtain direct structural evidence of the hexatic phase via scattering 
experiments. 

Effect of Multi-Layers on KTHNY Melting. The scaling of the order-parameter cor- 
relation functions in the simulations (figure 2) are consistent with the KTHNY behavior, 
implying that it is the vortex excitations that govern the equilibrium behavior in the quasi- 
two-dimensional systems and that the melting transition is defect-mediated. Moreover, for 
a quasi- two-dimensional monolayer, the Kosterlitz-Thouless transitions are continuous 37 , 
while for quasi-two-dimensional bilayers the Kosterlitz-Thouless transitions become first- 
order (figure 3). We ascribe this deviation from 2-d behavior to the interactions between 
the defect configurations in different layers, based on the following scaling arguments. 

On heuristic grounds 3 , the free energy of exciting a vortex pair of opposite wind- 
ing numbers relative to the ordered phase in the xy model in 2-d is given by 
AF = (irJ — 2k B T) log(L), J being the interaction energy of two neighboring spins with 
the same alignment and L being the distance between the vortex cores. This is the 
essence of the Kosterlitz-Thouless (KT) transition. If we consider two planes of xy mod- 
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els interacting with each other, and the vortices in one layer are perfectly in alignment 
with those in the second layer, the free energy relative to the ordered phase is given by 
Aj p(KT) = ( 27r j _ 2k B T) log(L), which is qualitatively the same as the KT behavior. This 
configuration, however, corresponds to a reduced entropy situation because the number of 
different ways of placing the aligned vortex pairs is the same as the single layer case. A 
second higher entropy scenario exists, where a vortex in one layer is aligned with a vortex 
of the opposite winding number in the other layer, with the cores displaced by distance A. 



Such a situation corresponds to a free energy (see Appendix): Appendix 

AF = AF KT + (tt - 1) J'L 2 + J'AL In L - k B T ln(l + nA 2 ) (21) 

where J is the interaction of nearest neighbors of spins with the same alignment, L is the 
in-plane distance between vortex cores, AF^ KT ) = {2nJ — 2k B T) log(L) is the Kosterlitz- 
Thouless free energy 3 , J' is the nearest-neighbor interaction between layers and in general 
can be different from J, and A is the "offset-distance" between vortices as defined above. 
Treating L as the order parameter, the free energy profile is obtained by the locus of points 
that minimize AF, the variational parameter being A. The locus of points minimizing AF 
is qualitatively different from that for the true KT behavior, given by AF^ KT ) (figure 9) 
. In particular, at small separations of the vortices, the bound state actually exists as Fig. 9 
a metastable state for T > T c , a clear signature of a first-order phase transition. Figure 9 
depicts the free energy profile for the monolayer case and bilayer case at a temperature k B T = 
1.5, assuming 2n J = 1 and (tt — 1) J' = 0.2 (where J, J' are in units of k B T). The plot marked 
by the solid line represents AF^ KT \ the true KT scenario, at a temperature, T, greater than 
the transition temperature, T c , for which the lowest free energy state corresponds to L = oo, 
i.e., the vortices are unbound. Moreover, the bound state (L = 1) is unstable and remains 
so, as long as T > T c , and therefore the transition is continuous (second-order). The set of 
dotted lines are plotted according to equation (21) for different values of the offset distance 
A. It is immediately clear that the locus of points minimizing AF is qualitatively different 
from the KT plot (equation for AF^ KT )). In particular, at small separations of the vortices, 
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the bound state actually exists as a metastable state for T > T c , a clear sign of a first- 
order phase transition. In other words, the vortex excitations that lead to equation A- 7 (see 

Appendix) imply a pathway for nucleation of the ordered phase, i.e., by the annihilation 

— » — * 

of vortices U + and U~ across layers. This is possible only if hopping of particles between 
layers are allowed, such as in our confined fluid system. 

The simulations and above arguments pertain to smooth pore walls. Therefore, the effect 
of the six-fold substrate potential was ignored while taking into account the strong potential 
due to the porous matrix. To our knowledge, this is the first report of the effect of multi- 
layers on KTHNY behavior. 



Effect of Strength of the Pore Potential. It was remarked earlier, based on the evidence 
from simulations and experiments, that the hexatic phase was stable over a large temperature 
range for CC1 4 and aniline. The ratio of the strength of the fluid-wall interaction to the fluid- 
fluid interaction, a (equation 3), has a significant effect on the melting phase diagram 33 (see 



figure 10); in particular, the larger the value of a, the larger the temperature range over Fig. 10 
which the hexatic phase is stable. Therefore, fluids with purely dispersive interactions 
(e.g., CCI4, methane, benzene etc.) confined in ACF, with large values of a show a large 
temperature range of stability for the hexatic phase. For strongly dipolar fluids, a is smaller, 
the temperature range of stability of the hexatic phase decreases 33 , and may vanish for 
sufficiently small a. 
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under grant number 2 P03B 01424. The experiments were performed in Adam Mickiewicz 
University, Poznan, Poland. The supercomputing time was provided by the San Diego Super- 
computing Center under a NRAC grant (number MCA93S011P). International cooperation 
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APPENDIX 



We give here the derivation of equation 21 for a bilayer of the xy model of spins on 
a lattice. We consider two planes of xy models interacting with each other; vortices with 
opposite winding numbers in one layer are perfectly in alignment with those in the second 
layer (i.e, the vortex cores span two layers). The free energy relative to the ordered phase 
is given by 3 : 

af (kt) = (2ttJ _ 2k B T) log(L) (A-l) 

which is qualitatively the same as the KT behavior (solid line in figure 9). This configuration, 
however, corresponds to a reduced entropy situation because the number of different ways 
of placing the aligned vortex pairs is the same as the single layer case. A second, higher 
entropy scenario exists, where a vortex in one layer is aligned with a vortex of the opposite 
winding number in the other layer, with the cores displaced by distance A. 

Vortices of positive and negative winding numbers whose cores are displaced by length A 

— * — * 

are represented as vector fields U + and U~ , given by: 

U+ = ( x - A)e x - ye y (A-2) 
U~ = -ye x + (x - A)e y (A-3) 

where e x , e y are unit vectors in the x and y directions. The interaction energy E relative to 

— * — * 

the ordered phase, for the given spin fields U + and U~ in each layer, is given by: 

E = J' [ [ U + -TJ- dxdy +nJ'L 2 (A-4) 
Jo Jo 

= (tt - l)J'L 2 + J'AL In L (A-5) 

and the excess entropy associated with placing the cores with an offset distance A is given 
by S = k B ln(l + nA 2 ). The free energy relative to the ordered state is given by: 

F = 2nJ\nL — 2k B T In L + (tt - 1) J'L 2 + J' AL\n L - k B T \n(l + n A 2 ) (A-6) 
= AF KT + (tt - I) J'L 2 + J'AL In L - k B T\n{\ + nA 2 ) (A-7) 
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FIGURES 

FIG. 1. Orientational correlation functions for the confined bilayer of adsorbate (H = 1.41 nm): 
liquid (360 K), hexatic (335 K), and crystal (290 K) phases. 

FIG. 2. Scaling behavior of orientational correlation functions. The solid line corresponds to 
zero slope, the dashed line corresponds to a slope of —1/4 and the dotted line corresponds to a 
slope of —2. 

FIG. 3. (a) Landau free energy functions for two different system sizes 60crff and 180<7ff for a 
pore width of H = 1.41 nm at the hexatic-crystal transition. The temperature was T = 293 K for 
the 60<7ff system, and T = 290 K for the 180<Tff system, (b) Corresponding results for H = 0.91 nm 
system at the liquid- hexatic transition, for the 60a // system (T = 390 K) and I8O07/ (T = 387 K). 

FIG. 4. The grand free energy per mole for liquid, hexatic and crystal phases {H = 1.41 nm); 
the cross-over points correspond to phase transitions and give the transition temperatures (see 
Table II). 

FIG. 5. DSC scan for CCI4 and aniline confined in activated carbon fiber ACF A-10 at a 
temperature scanning rate of 0.1 K/min. The peaks are interpreted as the liquid- hexatic (high 
T) and hexatic-crystal (low T) transitions. The approximate transition temperature is located at 
the tail of each peak in the low temperature side (see Figure 9.4.3 on page 550 of Chaikin and 
Lubensky 49 ). 

FIG. 6. DRS temperature scan (the sample was equilibrated at each temperature) for aniline 
confined in ACF-A10, as a function of temperature indicating phase transitions (dashed lines), (a) 
capacitance; (b) molecular orientational relaxation time (r). The disappearance of the relaxation 
branch due to conductance of the sample is believed to correspond to the liquid-hexatic transition, 
and the change in the relaxation time of the confined fluid is interpreted as the occurrence of the 
hexatic-crystal phase transition 52 . The relaxation branch corresponding to the Maxwell- Wagner 
effect occurs as a result of signal dispersion due to the suspended ACF particles in bulk aniline 52 . 
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FIG. 7. NDE for (a) CCU, and (b) aniline, confined in activated carbon fiber ACF A-IO. The 
peaks correspond to the liquid-hexatic (L/H) and hexatic-crystal (H/C) transitions. Lines are 
drawn through the points as guide to the eye. 

FIG. 8. NDE scaling for aniline confined in activated carbon fiber ACF A-10. The filled circles 
correspond to liquid above the L/H phase transition, and the filled squares correspond to the 
hexatic above the H/C phase transition. The solid lines represent a fit to the scaling law with 
the theoretically predicted exponents (see legend). For the liquid (filled circles) A = exp(34.52), 
B = 6.48, v = 0.5, T c = 42°C; for the hexatic (filled squares) A = exp(32.46), B = 0.338 , v = 0.37, 
T c = 27°C. 

FIG. 9. Free energy scaling for a set of interacting vortex pairs in a bilayer according to the 
equation for AF^^ and equation 21. The distances L and A are in dimensionless units (scaled 
by the lattice length) and the free energy is in units of ksT. 

FIG. 10. Global phase diagram of a fluid in a slit pore of width H = 3<7g- from simulations 
and free energy calculations 33 ' 37 (open symbols) and experiment (filled symbols) 33 ' 37 ' 55 ' 56,52 ' 41 ' 57 . 
Three different phases are observed: liquid (L), hexatic (H), and crystalline (C). The dashed 
line represents an extrapolation of the phase boundaries based on MC simulations without free 
energy calculations. The simulations are for a LJ fluid in slit-pore with different values of a. The 
experiments are for various adsorbates confined within activated carbon fibers (ACF, mean pore 
width 1.4 nm): H 2 55 {a = 0.51), C 6 H 5 N0 2 56 (a = 1.22), C 6 H 5 NH 2 52 (a = 1.75), CH3OH 55 
(a = 1.82), CCI4 41 ' 37 (a = 1.92), C 6 H 6 57 (a = 2.18). 
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Table I: Ginzburg parameter 
for different L and T. 



T/K 10CT ff 40a ff 60a s 180a ff 
340 8.98 1.95 0.117 0.038 
330 6.4 1.7 0.102 0.029 
320 4.29 1.17 0.098 0.019 
310 3.11 0.97 0.16 0.0189 
300 3 0.937 0.176 0.0189 



Table II. Transition temperatures from simulation and 



experiment for CCU and aniline in ACF-10. 



Fluid 


Low T c /K (H/C) 


High T c /K (L/H) 


Simulation 


Expt. 


Simulation 


Expt. 


CCI4 
Aniline 


290 


298 a , 295 6 - c 
301 a , 298 6 , 300 c 


348 K 


348 a , 353 6 , 352 c 
325 a , 324 6 , 315 c 
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